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This study demonstrates that coupling of a material thermal response code and a flow solver using 
finite-rate gas/surface interaction model provides time-accurate solutions for multidimensional ablation 
of carbon based charring ablators. The material thermal response code used in this study is the Two- 
dimensional Implicit Thermal Response and Ablation Program (TITAN), which predicts charring material 
thermal response and shape change on hypersonic space vehicles. Its governing equations include total 
energy balance, pyrolysis gas momentum conservation, and a three-component decomposition model. 
The flow code solves the reacting Navier-Stokes equations using Data Parallel Line Relaxation (DPLR) 
method. Loose coupling between material response and flow codes is performed by solving the surface 
mass balance in DPLR and the surface energy balance in TITAN. Thus, the material surface recession is 
predicted by finite-rate gas/surface interaction boundary conditions implemented in DPLR, and the 
surface temperature and pyrolysis gas injection rate are computed in TITAN. Two sets of gas/surface 
interaction chemistry between air and carbon surface developed by Park and Zhluktov, respectively, are 
studied. Coupled fluid-material response analyses of stagnation tests conducted in NASA Ames Research 
Center arc-jet facilities are considered. The ablating material used in these arc-jet tests was a Phenolic 
Impregnated Carbon Ablator (PICA). Computational predictions of in-depth material thermal response 
and surface recession are compared with the experimental measurements for stagnation cold wall heat 
flux ranging from 107 to 1100 W/cm 2 . 


Nomenclature 

= mass fraction 
= diffusion coefficient, m 2 /s 
= enthalpy, J/kg 

= thermal conductivity of translational temperature, W/m-K 
= thermal conductivity of vibrational temperature, W/m-K 
= mass flux, kg/m 2 -s 
= molecular weight, kg/mole 
p = pressure, N/m 2 

q com , = convective heat flux at surface, W/m 2 

q cw = In-depth conductive heat flux at surface, W/m 2 

q nv = radiative heat flux at surface, W/m 2 

r c = corner radius, cm 

r n = nose radius, cm 

^max = maximum allowed change on surface recession, m 
T = temperature, K 

T 

max = maximum allowed change on surface temperature, K 
T t = translational temperature, K 

T v = vibrational temperature, K 
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T 

J 00 

= environment temperature, K 

V 

= mass injection velocity, m/s 

w 

= source term of gas-surface interactions, mole/m 2 -s 

a 

= absorptance 

£ 

= emissivity 

P 

= total density, kg/m 3 

a 

= Stefan-Boltzmann constant, W/m 2 -K 4 

subscripts 

c 

= char 

g 

= pyrolysis gas 

i 

= gas species 

w 

= wall 


I. Introduction 

Vehicles designed for Earth entry at superorbital velocities, as well as those designed for ballistic entry at orbital 
velocities, typically use thermal protection system (TPS) materials that pyrolyze and ablate at high temperature for 
mass-efficient rejection of the aerothermal heat load. For design and sizing of ablating TPS materials, it is 
imperative to have reliable numerical procedures which can accurately compute both the aerothermal environments 
and the material surface ablation and internal thermal response. It has been demonstrated that accurate prediction of 
ablative heat flux also requires fliuid-solid shape change coupling simulation. 1 In our previous fluid-solid shape 
change coupling computations, the ablative carbon surface was assumed to be at chemical equilibrium. 2 Chemical 
equilibrium assumption is good for many space entry applications, but may not be true for all conditions. 
Additionally, the aerothermal environments computed by Data Parallel Line Relaxation method (DPLR) 3 were those 
for a fully catalytic non-ablating surface. Thus, an engineering correlation with blowing reduction parameter had to 
be introduced in the Two-dimensional Implicit Thermal Response and Ablation program (TITAN) 4 simulation to 
take into account the effect of mass injection on reduction of convective heat flux. A mass transfer coefficient also 
had to be defined based on the heat transfer coefficient under the assumption of relatively thin boundary layer for the 
computation of char recession rate. 

The purpose of this paper is to present a new TITAN-DPLR coupling simulation system and its applications to 
carbon based charring ablators. The boundary condition with general finite-rate chemistry for gas/surface 
interactions recently implemented in the DPLR code by MacLean 3 is used to predict char mass injection rate. Thus, 
surface thermal chemistry coupling is also introduced in this simulation system, in addition to surface shape change 
coupling. In this simulation system, the surface species mass balance is performed in DPLR, and the surface energy 
balance is performed in TITAN. The chemical equilibrium assumption can then be removed. The hot-wall ablating 
convective heat flux is directly computed in DPLR based on the surface temperature and pyrolysis gas injection rate 
computed in TITAN. The blowing reduction correlation is no longer needed in the material response simulation. 
Two sets of finite-rate gas/surface interaction chemistry between air and carbon surface are studied. They were 
developed by Park 6 and Zhluktov, 7 respectively. Coupled fluid-material response analyses of stagnation tests 
conducted in NASA Ames Research Center arc -jet facilities are performed. The ablating material used in these arc- 
jet tests was a Phenolic Impregnated Carbon Ablator (PICA). 8 Computational predictions of in-depth material 
thermal response and surface recession are compared with the experimental measurements for stagnation cold wall 
heat fluxes from 100 to 1100 w/cm 2 . 

II. Surface Thermal Chemistry Coupling 

The interface between a TPS material and its surrounding flow field can be defined by solving species mass 
conservation and energy balance equations. Species mass conservation at the surface of TPS material is written as: 9 

-pD l VC l + pv w C t = M, w, + m g C ig ( 1 ) 

The first term on the left-hand side is mass transfer through diffusion, and the second term is mass transfer due to 
convection. On the right-hand side are the source terms due to gas-surface interaction and pyrolysis gas injection. 
Based on global mass balance at the surface, the following equation for the total mass blowing rate is expressed: 
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( 2 ) 


pv w = m c + m g 

The total convective heat flux to the surface for the flow field with two-temperature model is given as: 

Qconv = ~ k t yT , ~ W + Yj h iP D i SJC i • ( 3 * ) 

Energy conservation equation at the surface is written as 

Qconv T rti c (Ji c hw') T h w ) T tx w q rw cj£ w (T w T^) q cw 0 (4) 

The first term in Equation (4) is the total convective heat flux, the second and third terms represent the heat of 
ablation, the fourth and fifth terms are radiation absorption and emission, respectively, and the final term is the rate 
of heat conduction into the TPS material. 

To obtain the solutions for Equations (1) to (4) requires computations of thermal and species diffusion rates of 
flow field at the surface as well as thermal diffusion and pyrolysis gas injection rates of TPS material at the surface. 
In a strongly coupled simulation, the governing equations for both fluid and solid are solved simultaneously along 
with Equations (1) through (4). Elowever, for many simulations, flow solver and material response code are two 
independent programs, and a loosely coupled simulation needs to be applied. In a loosely coupled simulation, the 
governing equations for fluid and solid are solved separately, and Equations (1) to (4) have to be solved either in the 
flow code or in the material response code. Thus, communication between two codes needs to be established for 
exchanging information on surface thermal chemistry and shape change. If the material surface is at chemical 
equilibrium, both the mass and energy conservation equations at the surface can be performed in the material 
simulation code as described in our previous work. For a chemical equilibrium surface, the chemical species at the 
surface are determined by using a chemical equilibrium code, such as ACE 10 and MAT. 11 In our chemical 
equilibrium analyses, the hot wall ablating heat flux was estimated from the cold wall non-ablating heat flux using 
an engineering correlation with a blowing reduction parameter. Also, species mass transfer rate was assumed to be 
proportional to heat transfer rate based on a constant Lewis number. Under these assumptions, the surface thermal 
chemistry coupling is not required, and only shape change information has to be exchanged between flow solver and 
material response code to correctly predict the aerothermal environments and material thermal response. For a 
general finite-rate surface boundary condition, the approach taken here is having the species mass conservation 
equation solved with flow-field governing equations and total energy balance equation solved with solid material 
governing equations. Surface thermal chemistry and shape change have to be shared between flow solver and 
material solver in a time lagged manner. The schematic diagram in Figure 1 depicts how this loosely coupled 
simulation is performed. This simulation must be time accurate and starts from the material response simulation 
code using the boundary conditions of cold wall non-ablating convective heat flux and pressure estimated based on 
initial free stream conditions. As the maximum surface temperature change or surface recession exceeds the pre- 
determined value, material response computation is temporarily put in a waiting mode and the latest predicted 
surface temperature and pyrolysis gas injection rates are used as the boundary conditions in performing a flow field 
simulation. The pyrolysis gas is assumed to be at chemical equilibrium before injected into the adjacent air. This 
assumption was proved to be reasonable for conditions studied in this paper. 12 The computational grid system for 
flow field is reconstructed based on the shape change predicted by the material response code, and the free stream 
conditions are also updated according to the flight trajectory or arc heater condition. Each flow simulation is a 
steady-state computation. When the steady-state flow solution is obtained, the time -dependent material response 
simulation is resumed using the newly predicted surface heat flux, pressure, enthalpy, and char mass injection rate as 
the boundary conditions. Again material response simulation continues until maximum surface temperature change 
or surface recession exceeds the maximum allowed value. This process repeats until the end of flight trajectory or 
arc -jet exposure. 
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Figure 1. Surface thermochemistry and shape change coupling between FIAT/TITAN and DPLR. 


III. Chemistry Models 

Two sets of gas/surface interaction model between air and carbon surface are studied in this work. The first set 
was developed by Park 6 and the second set was developed by Zhluktov. 7 In this study, the nitridation reaction in the 
original Park’s model is replaced by the nitrogen recombination reaction to be consistent with the observation in the 
arc -jet tests using pure nitrogen gas. Zhluktov’s model is more sophisticated as compared with Park’s. In Zhluktov’s 
model, the numbers of empty surface sites, E(s), and sites occupied by oxygen, O(s), or nitrogen, N(s), are 
considered to predict both the forward and backward reaction rates. Park’s model uses an engineering approach 
considering only the forward reaction. Here C(b) represents the solid carbon species. The gas/ surface reactions used 
in this study are listed below, and their rates can be found in Refs 6, 7, and 9: 

Modified Park: 

#1 O + C(b) -> CO 
#2 0 2 + 2C(b) -» 2CO 
#3 2N -> N 2 
#4 3C(b) -» C 3 
#5 C 3 ^ 3C(b) 

Zhluktov: 

#1 O + E(s) <-> O(s) 

#2 20(s) <-* 0 2 + 2E(s) 

#3 0 2 + E(s) *-> O + O(s) 

#4 C0 2 + E(s) <-> CO + O(s) 

#5 O(sl) + C(b) <-> CO + E(sl) 

#6 O + O(s) + C(b) <-> C0 2 + E(s) 

#7 20(s) + C(b) <-> C0 2 + 2E(s) 

#8 C + E(s) ~ E(s) + C(b) 

#9 C 2 + 2E(s) ~ 2E(s) + 2C(b) 

#10 C 3 + 3E(s) <-> 3E(s) + 3C(b) 

#11 N + E(s) «-» N(s) 

#12 N 2 + E(s) + N(s) 

Marschall developed a general formulation for finite -rate gas and surface interactions. 13 The DPLR code was 
enhanced by MacLean based on Marschall’ s work to solve the surface species mass balance equation for the fmite- 
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rate gas surface reactions. The detail of his implementation can be found in Ref. 5. All flow simulations presented in 
this work are performed using this version of DPLR code. 

There are twenty two gas phase chemical species used in this study for the simulation of PICA and air- Argon 
interactions. None of the ion species are considered. The chemical species are CO 2 , CO, N 2 , 0 2 , NO, C 2 , C 3 , CN, H 2 , 
HCN, C, N, O, H, CH, CH 2 , C 2 H, C 2 H 2 , C 4 , C 3 H, C 4 H, and Ar. These species were selected based on the chemical 
equilibrium computation for a PIC A/air- Argon system. Their enthalpy makes up more than 95% of total gas 
enthalpy as compared with the baseline 119 chemical equilibrium species model developed by Orion TPS Advanced 
Development Project. 14 The gas phase chemical reactions considered in the simulations are 

#1 C0 2 + M CO + O + M 

#2 CO + M<->C + 0 + M 
#3 N 2 + M~N + N + M 
#4 0 2 + M~0 + 0 + M 

#5 N0 + M<->N + 0 + M 
#6 C 2 + M~C + C + M 
#7 C 3 + M<->C 2 + C + M 
#8 CN + M^C + N + M 
#9 H 2 + M~H + H + M 
#10 N 2 + Ar N + N + Ar 
#11 0 2 + Ar <->0 + 0 + Ar 
#12 NO + Ar <-> N + O + Ar 
#13 N0 + 0~0 2 + N 
#14 N 2 + 0~N0+N 
#15 C0 + 0~0 2 + C 
#16 C0 2 + 0~0 2 + C0 
#17 C0 + C~C 2 + 0 
#18 CO + N <-> CN + O 
#19 N 2 + C<->CN + N 
#20 CN + O NO + C 
#21 CN + C~C 2 + N 
#22 HCN + H~CN+H 2 
#23 CH + M^C + H + M 
#24 CH 2 + M<->C + H 2 + M 
#25 CH 2 + M«->CH + H + M 
#26 C 2 H + M~C 2 + H + M 
#27 C 2 H 2 + M<->C 2 H + H + M 
#28 C 2 + C 2 C 3 + C 
#29 C 2 + H 2 ~C 2 H + H 
#30 CH 2 + C CH + CH 
#31 CH 2 + CH 2 <-> C 2 H 2 + H 2 
#32 CH 2 + C C 2 H + H 
#33 CH 2 + C 2 H «-» CH + C 2 H 2 
#34 CH + CH <-> C 2 H + H 
#35 CH + C 2 H <-> C 2 H 2 + C 
#36 CH 2 + CH 2 «-» C 2 H 2 + H + H 
#37 C 2 + C 2 + M~C 4 + M 
#38 C + CH <-> C 2 + H 
#39 C + C 2 H~C 3 + H 
#40 C 2 + CH <-> C 3 + H 
#41 C 2 + C 2 H <-> C 4 + H 
#42 CH + H<->C + H 2 
#43 CH 2 + H <-> CH + H 2 
#44 C 2 H + C 2 H <-> C 2 H 2 + C 2 
#45 C + C 2 H 2 ~C 3 H + H 
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#46 C 2 + C 2 H 2 <-» C 4 H + H 
#47 C 2 + C 4 H <-> C 2 H + C4 
#48 C 2 H + C 2 H C 4 H + H 
#49 C 3 H + H~C 3 + H 2 
#50 C 4 H + H~C 4 + H 2 
#51 H + C 4 + M~C 2 H + M 
#52 CH + CH~C 2 + H + H 
#53 C 2 H + H 2 C 2 H 2 + H 

The reaction rates for reactions #1-22 are taken from the work of Olynick et al. 15 for Stardust earth entry 
simulation. The rates for reactions #23-36 are taken from Go keen’s 16 paper for simulations of Titan atmosphere 
entry. The rates for reactions #37-53 are taken from the study of Kruse et al. 17 for high-temperature pyrolysis of 
Acetylene. 


IV. Results 

The computations presented in this section focus on the analysis of surface recession and in-depth thermal 
response for the PICA heat-shield material. Model validation is accomplished by comparison of predictions with 
data from seven selected arc -jet tests conducted over a range of stagnation heat flux and pressures ranging from 107 
W/cm 2 at 2.3 kPa to 1100 W/cm 2 at 84 kPa. The test conditions for these seven test cases are listed in Table 1. 
Figure 2 shows the model shape used in these tests. The nose radius equals to the model diameter, and the sides are 
cylindrical. This is so called “iso-q” geometry (r„ = 5.08 cm, and r c /r n = 1/16). 


Table 1: Seven arc-jet cases selected for detailed analysis. 


Case number 

Stagnation Point heat flux, 
W/cm 2 

Stagnation point pressure, 
kPa 

Exposure time, 
s 

1 

107 

2.3 

55 

2 

169 

5.0 

60 

3 

246 

8.5 

42 

4 

395 

17.2 

34 

5 

552 

27.3 

30 

6 

744 

31.0 

27 

7 

1102 

84.4 

10 
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Figure 2. Model geometry and material map. 


Figures 3a and 3b present the predicted histories of mass injection rates, heat flux, and temperature at stagnation 
point for test case 4 using modified Park’s surface chemistry. Char mass injection rate and surface heat flux are 
predicted by DPLR (red lines). Pyrolysis gas mass injection rate and surface temperature are predicted by TITAN 
(blue lines). The distributions of mass injection rates, heat flux, and temperature along the model surface at the end 
of exposure time (t = 34 s) are shown in Figs. 4a and 4b. The predicted stagnation-point char mass injection rate is 
almost identical to data (symbol), which is the measured stagnation point recession multiplied by char density and 
divided by exposure time. Figure 4c shows major species mass fractions (CO, NO, C 3 , FF, N 2 , N, O, and FI) along 
stagnation stream line. The comparison of in-depth thermal response at various depths between computation (black 
lines) and TC data (red lines) is presented in Fig. 5. The thermocouple locations within arc-jet models are listed in 
Table 2. The TC locations for case 4 follow Option B. The agreement between TC data and computation is very 
good. 

The comparison of stagnation point heat flux computed based on different surface chemistry models is shown in 
Fig. 6a. The predictions using modified Park’s model (squares) agree well with the chemical equilibrium surface 
model with blowing reduction parameter equal to 0.5 (circles). The heat fluxes predicted using Zhluktov’s surface 
chemistry (triangles) appear to be lower than those using Park’s model at high heating conditions (cases 4-7). Figure 
6b is the comparison of stagnation-point char mass injection rates. Generally speaking, char recession rates 
computed using Modified Park’s model are in excellent agreement with data except case 7, in which the prediction 
is about 20% below the measurement. Zhluktov’s surface chemistry generally under-predicts char recession rates for 
all the seven cases. 

In our proposed paper, the computations using both Park’s and Zhluktov’s air-carbon interaction models for all 
seven cases will be presented. The detailed comparisons between current finite-rate approach and the chemical 
equilibrium approach adopted in our previous paper 3 will also be discussed. 


Table 2: Thermocouple locations within arc-jet models. 



Axial depth, 

cm 


TC number 

Option A 

Option B 

Option C 

1 

0.381 

0.508 

0.635 

2 

0.762 

0.889 

1.016 

3 

1.143 

1.270 

1.397 

4 

1.524 

1.651 

1.778 

5 

3.048 

3.048 

3.048 
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Stag Point Temperature (C) = Stag Point Mass Injection Rate (kg/m 2 -s) 



3a: Stagnation point mass injection rates for case 4 using Park’s model. 



450 

400 

350 ^ 
<_> 

300 W 

X 

250 

4 — ' 

CD 

200 

+-» 

150 '§ 

O- 

0 £> 

100 is 

i/i 

50 

0 


Figure 3b: Stagnation point heat flux and temperature for case 4 using Park’s model. 
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Figure 4a: Surface temperature and heat flux distributions using Park’s model for case 



Figure 4b: Mass injection rate distributions using Park’s model for case 4. 
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Temperature (C) 



Figure 4c: Species mass fractions along stagnation stream line for case 4. 
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Figure 5: Comparison between prediction and TC data for case 4 using Park’s model. 
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Figure 6a: Stagnation point heat fluxes. 
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Figure 6b: Stagnation point char mass blowing rates. 
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